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Abstract 

This paper presents how the space of spheres and sheUing may be 
used to delete a point from a d-dimensional triangulation efficiently. In 
dimension two, if k is the degree of the deleted vertex, the complexity 
is O(felogfc), but we notice that this number only applies to low cost 
operations, while time consuming computations are only done a linear 
number of times. 

This algorithm may be viewed as a variation of Heller's algorithm 



|Hel90, Mid93|, which is popular in the geographic information system 
community. Unfortunately, Heller algorithm is false, as explained in this 
paper. 



1 Introduction 

The computation of the Delaunay triangulation of a set S of n points in 
the plane is one of the classical problems of computational geometry. 

Many structures and algorithms have been proposed in the past to 
compute Delaunay triangulations. Some of these algorithms have the two 
following properties: they are incremental and they do not used compli- 
cated data structures in addition to the triangulation itself. Among these 
algorithms, let us mention the historica l algori t hm of Green and Sibso n 
|GS78[, or some other variants ]MSZ96|, [BD95[, pev98|, [DLM9E|, |Lem97l. 



All perform a walk in the triangulation to accelerate point location. 

The advantage of that category of incremental Delaunay algorithms is 
that they may easily be turned into fully dynamic Delaunay algorithms. 
Since there is no complicated data structure for point location, the deletion 
of a point is reduced to the deletion in the triangulation itself. 

Definition and notations 
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Given a set S of points in d- dimensional space, T>T{S), the Delaunay 
triangulation of S is defined by the following property: d + 1 points of S 
are the vertices of a Delaunay simplex if and only if the sphere passing 
through these points does not contain another point of S in its interior 
(see Figure |l| for a Delaunay triangulation in two dimensions). 

Given the Delaunay triangulation VT{S) and a vertex p in VT{S), 
we address the problem of finding T>T{S \ {p})- 

In two dimensions, the natural parameter to evaluate the complexity 
of this problem is the degree fc of p in T>T{S), since the deletion of p means 
that k triangles must be removed from the triangulation and k — 2 new 
triangles must be created to fill this hole. In the worst case, k may be \S\, 
but if p is chosen randomly in S, then it is well known that the expected 
value of k is 6, without any assumption on the point distribution. 

In higher dimensions, the number of simplices incident to p is not 
directly related to the number of simplices created to fill the hole. We 
will let / denote the sum of these two numbers. In the worst case, the 
whole Delaunay triangulation may be affected, and / = 0{n^ 2 J). This 
distribution does not reflect practical conflgurations, and a constant value 
of / is more likely in practice. For a uniform point distribution, the 
expected value of / may be shown to be constant. In three dimensions, 
the expect ed numb er of deleted tetrahedron for Poisson distribution is 
fyr^ ~ 27 |OBS92 |. 

Thus, even if we do not want to neglect the possible case of a big value 
for k, we have to keep in mind that a good algorithm must perform well 
on small values of k. 

Previous related work 

Classical Computational Geometry has already addressed the problem 
of deleting points from Delaunay triangulations. This can be done wit h 
optimal asymptotic complexity 0{k) in 2 dimensions [Che86, AGSS8£| ]. 
But these algorithms are a little bit too intricate and the big O hides too 
important a constant for them to be good algorithms for reasonable values 
of k (we give this constant in Section 4.1 for Chew's algorithm, Aggarwal 
et al's algorithm is even more complicated). 

Practitioners often prefer algorithmic simplicity to theoretical optimal- 
ity, and favour a simple suboptimal O(fc^), implementation of the deletion 
algorithm. This may be achieved, for example, by flipping, to reduce the 
degree of the deleted vertex to 3, and flipping again to restore the Delau- 
nay property. Another simple algorithm consists in finding the Delaunay 
triangle incident to an edge of the hole in 0{k) time, which also yields an 
O(fc^) time algorithm. 

A very simple O(felogfe) solution was suggested by Heller [Hel90, 



Mid93|, in which successive ears are filled in turn. In that way, dur- 



ing the algorithm we always have a simple polygon of decreasing size to 
triangulate. Unfortunately, this solution is wrong, but we will show in 
this paper how to correct it. 
Overview 

In this paper, we provide a very simple and efficient O(fclogfc) al- 
gorithm to delete a v ertex in a pla nar Delaunay tr iangulation based on 
shelhng pM7l\ ^ei86| and duality |DMT92| , |Ped70| . We also discuss the 
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effective complexity of this algorithm and of a few others for small values 
of k. We will study the different kinds of geometric predicates necessary 
for these algorithms. 

This algorithm generalizes well in higher dimensions: its time com- 
plexity becomes 0{f log f), where / is the number of tetrahedra created, 
and it also generalizes to regular triangulations (power diagrams). 

2 Two dimensional algorithm 

A deletion algorithm has to remove all triangles incident to p and retrian- 
gulate "Delaunay-wise" the star-shaped polygon H = {go, qi, . . . , qt-i, Qk = 
go} created by these removals (see Figure Ol). 

2.1 Ears, and a wrong algorithm 

We first define what an ear of a polygon is. Three consecutive vertices 
qiqi+iqi+2 along //'s boundary are said to form an ear of H if the line 
segment qiqi+2 is inside H and does not cross its boundary. An ear of 
H is said to be Delaunay, if the circle through qi gi+i and qi+2 does not 
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Figure 2: The smallest potential ear may not belong to the final Delaunay 
triangulation. The shaded triangle is the ear with smallest circumradius, but 
its associated circle contains (73 and thus is not a Delaunay ear. 



contain an y other vertices of H in its interior. Heller [ |Hel90| (also cited 



by Midtb0 [ Mid93[ |) claimed (without proof) that among all the potential 



ears qiqi+iqi+2 of H, the one having the circumcircle with smallest radius 
is a Delaunay ear. This claim is false, as illustrated by Figure ^: on the 
left handside are shown the Delaunay triangulation 'DT{S) and the hole 
H to be retriangulated; on the right handside are shown two potential ears 
qoq\q2 and ^15293. qoqiq2 has the smallest circumcircle among all ears of 
H, but it contains qs, which invalidates Heller's claim. Heller's mistake 
is to assume that when we deform a circle through §152 to maximise its 
portion inside H, the radius increases, but this is true only if the center 
of the circle is inside H, which is not the case in Figure ^. 

In fact, the idea of finding an ear belonging to the Delaunay triangu- 
lation works, but with another criterion, as explained below. 

2.2 Delaunay and convex hull 

There exists a well known duality between Del aunay t r iangulat ions in 
dimension d and convex hulls in dimension d + 1 lAur87| , |pMT92|l . If we 



associate to a point p = {x,y) £ S a point p* = {x,y,x + y ) on the 
paraboloid H of equation z = + y^, the Delaunay triangulation of S 



4 



is the projection of the convex hull of the 3D points. The reason is that 
for p, q,r, s £ S , p is inside the circle Cqrs through qrs if and only if p* is 
below the plane Pqrs through q*r*s* (II n Pqrs projects onto circle Cqra)- 
We even have the equality between the power of p with respect to Cqrs 
and the signed vertical distance between p* and Pgrs-fj 



2.3 Shelling 



Convex hulls m ay be computed by the shelling algorithm [seiSCj. The 



shelling [ BM71 of a convex polyhedron P is the enumeration of its faces 
in some appropriate order. Imagine an observer is moving along a line 
I going through the polyhedron, starting at the intersection of P and I. 
At the starting position, the observer can only see one face of P (the 
face intersecting its trajectory) and when she moves away from P she 
discovers other faces one by one; at infinity, the observer sees "half" of 
P. Then the observer turns round at infinity on the opposite side of I 
(where she sees the other half of P), and then moves on I, enumerating 
the faces of P when they disappear from her view. This order is called 
the shelling order of P with respect to /; it has the property that the set 
of enumerated faces remains simply connected during the enumeration. 
Seidel's algorithm for convex hull reports the faces of the convex hulls in 
that order, by maintaining a priority queue of potential new faces. These 
new faces are of two kinds depending on whether the next face contains 
a new vertex, or else connects already visible vertices. 



2.4 Deletion in Delaunay 

Let us now use the idea of finding Delaunay ears to retriangulate the hole 
created by the deletion of a point p in T)T{S). Using duality with convex 
hulls, the problem is transformed into filling the hole in the convex hull 
created by the deletion of p*. Then we may note that if we use shelling 
order with respect to the vertical line through p* , an observer (going up) 
reaching p* sees the boundary of this hole exactly; thus, the end of the 
shelling procedure corresponds exactly to the triangulation of the hole. 
This partial shelling is easier to implement than Seidel's algorithm, since 
all the vertices are already visible and thus only one kind of potential new 
faces has to be found. The already noted correspondence between vertical 
distance and power yields the following lemma. 

Lemma: Consider a polygon H — {go, qi, ■ ■ ■ ,qk~i,qk ~ qo} 
and a point p such that the edges ofqiqi+i belongs to the Delau- 
nay triangulation of{qo,qi, qk-i,p}- If \power{p, circle{qi, qi+i, qi+2))\ 
is minimal, then qiqi+2 is an edge of the Delaunay triangulation 

of{qo,qi, ■ ■ ■,qk-i} 

^ If C is a circle of center x and radius r, p is a point and I is a line through p intersecting 
C in t and u, then power{p,C) = \xp\'^ — = ptpu where yz is the signed length of yz. 
power{p,C) is zero on C boundary, negative inside and positive outside. When C is given by 
three points, the power may be known without computing the circumradius, as explained in 
Section |3.l|. Power is negative inside the circle and positive outside. 
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Thus the deletion of a point may be implemented in a simple way, by 
maintaining a structure to store the ears. The ears are naturally ordered 
along the boundary of the hole, each with its priority. This structure must 
support the following operations : find next and previous ear according 
to counterclockwise order along H's boundary, delete ear with minimum 
priority and update the priority of an ear. This structure may be imple- 
mented with any dictionary structure augmented by next and previous 
pointers. 

Algorithm Delete (T>T(S),p)) 



1. Let qoqi ■ ■ ■ Qk-i be the vertices incident to p in 'DT{S) in ccw order 
around p; 

2. Let Q be a priority queue; 

3. for i = to fc — 1 

4. do ear ^qiqi+iqi+2; 

5. if counterclockwise(gi9i+iq'i+2) 

6. then p <— 00; //not an ear 

7. else p« power (p, ear); / /inside circle powerK Q 

8. Q. insert (p, ear); //insert(priority,key) 

9. while (5-size()> 3 

10. do ear <— Q.minimum(); 

11. create triangle ear and link it to its two existing neighbors; 

12. earO ^ear. previous; 

13. earl <— ear .next; 

14. ear0.vertex(2) <— ear. vertex(2); earO.next <— earl; 

15. earl.vertex(O) <— ear. vertex(O); earl.previous <— ear2; 

16. (5.delete(eor); 

17. (5.modify-priority(earO); 

18. (5.modify-priority(earl); 

19. ear <— Q.minimum(); //the three last ears are identical 

20. create triangle ear and link it to its three existing neighbors; 



Higher dimensions The generalization to d dimensions is easy. The 
boundary of the region to retriangulate is a simple polyhedron H, and 
the ears are simplices formed by the vertices of two incident facets of H. 
The difference is that the same simplex may correspond to 0{d^) pairs of 
incident facets, and that the creation of an ear may modify 0{d^) other 
ears in the priority queue. 

3 2D Analysis 
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3.1 Power computation 



An analytical expression of the power of p with respect to qoqiq^ is 
power (p, cirde{qo,qi,q2)) 
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The 3x3 determinant is the orientation test of qcqiq^z, and the 4x4 de- 
terminant the incircle test of p with respect to go5i52- First note that if 
qoqiq2 has the wrong orientation, then it is not an ear and thus the 4x4 
determinant does not need to be computed, and also that the orientation 
test is a minor of the incircle test, and thus the power computation just 
requires one division in addition to the usual incircle test. 

Using a dynamic programming development of the determinant, the 
power computation requires 14 additions, 15 multiplications and one di- 
vision. 

Finally, note that if only q2 changes, we do not need to recompute 
everything. In fact, we can do it with 6 additions, 7 multiplications and 
one division. 

3.2 Complexity 

The theoretical asymptotic complexity of this algorithm is clearly 0{k log k) , 
but this complexity only concerns the management of the priority queue, 
which involves relatively cheap operations (pointer manipulations and 
comparisons of computed powers). 

Since the most expensive geometric operations are the power compu- 
tations, we shall count the number of such operations exactly. The initial 
size of the priority queue is fc, and thus its initialization requires at most 
k power computations. Each ear creation implies the modification of two 
other ears, and thus two powers must be recomputed, and the deletion is 
completed when the size of the queue is 3; thus the total number of power 
computations is A; -|- 2{k — 4) = 3fc — 8. 

It is possible, as noticed above, to update the power of p with respect 
to qiqi+iqi+2 when ear qi+iqi+2qi+3 is processed. In the new polygon H\ 
{qi+2}, qiqi+iqi+3 is an ear and the power of p with respect to qiqi+\qi+3 
may be obtained by updating the power of p with respect to qiqi+iqi+2 at 
a cheaper cost. Then the total number of computations becomes 2k — 4 
power computations and fc — 4 power updates. 

3.3 Robustness issues and degeneracies 

The algorithm presented above does not address robustness issues. If 
floating point arithmetic is used to perform power computations, the re- 
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suits are rounded and their comparisons could be evaluated erroneously. 
We can first observe that the deletion algorithm will terminate even with 
incorrect arithmetic: it fill ears in turn and thus constructs a topological 
triangulation, which may be non Delaunay, or even have a non-planar 
embedding. But, even if producing a non-exact Delaunay triangulation 
may be acceptable for the deletion algorithm, it is unacceptable for many 
insertion algorithms which are not capable of processing non-exact trian- 
gulations. 

The usual way to solve robustness issues consists in using exact arith- 
metic. To ensure good performance, we can use arithmetic filters to use 
exact computations in power comparisons only in difficult cases where 
the two powers are close. Exact computations are then used to take the 
right decision. This approach of filtering out easy cases h as been proven 
efficient on the orientation and incircle tests IPPOSj |DP99 L Degeneracies 
can be solved using perturbation techniques [pei98^ |ADS97 . 



4 Alternative methods and practical re- 
sults 

4.1 Alternative methods in two dimensions 

Diagonal flipping One of the interests of a method using diagonal 
flipping is that it does not introduce new geometric predicates: it only uses 
incircle tests, and thus has lower degree and generalizes easily to various 
metrics. However such a method may require O(fc^) incircle tests in the 
worst case, and is not so simp le to code efficiently. A good implementation 



of a flipping method |Sno98| will retriangulate H by basically linking all 



qi to go and flipping the edges turning around go- In the worst case, such 
a method may use (A: — 3) -I- (A: — 4) -I- ... -I- 1 = (k-2)(k-3) jjj^jj-gjg ^gg^s. 

As this flipping method is not so easy to code, a simpler solution may 
consist in maintaining a queue of edges to be tested for potential flip. 
Each time the diagonal of a quadrilateral is flipped, the four edges of the 
quadrilateral are inserted in the queue. This simpler algorithm will make 
much more incircle tests, since the flip are not performed in some relevant 
order as above. 



Edge completion A second method consists in finding qi such that 
qoQiQi is a Delaunay triangle, which may be done in — 3 incircle tests, 
and triangulating recursively the two holes. In the worst case, the number 
of incircle tests is exactly the same that in the flipping method. 



Randomized algorithm In Chew |Che86] randomized algorithm, 
each point to reinsert requires an expected number of 5 incircle tests, 
which yields a total of 5k + 0{1). 



On alternative methods For small value of < 9, flipping or edge 
completion requires less incircle tests computations, but the simplicity of 
our "ear-queue" algorithm and its good performance for k > 9 make it a 
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very good candidate for Delaunay vertex deletion. Nevertheless, it can be 
interesting to treat as special cases some small k values such as fc = 4 and 
fc = 5. 



4.2 Higher dimensions 

In higher dimensions, flipping, edge completion and shelling algorithms 
generalize but things became more difficult. The flipping must be done i n 
a higher dimension, which makes it more intricate to implement | |ES9e| ]. 
Edge completion transforms into facet completion, and must deal with the 
triangulation of non simply connected polyhedra. Shelling is the easiest 
method to generalize. Furthermore, the increase in the average value of 
k with the dimension reinforces its advantage over alternative candidates 
in higher dimensions. 



4.3 Experimental results 

Code and data This algorithm was implemented within the author'; 



simple hierarchical structure [Dev' 



Robustness issues are solved using 24 bits integers to store points coor- 
dinates. Geometric predicates are computed with approximate arithmetic 
and the exactness of the result can be ensured by static and semi-static 
fihers^. The filters failure are backed up by exact computations. 

According to paragraphs above, we have coded several versions of the 
deletion procedure. 

• earS. The basic version using the ear queue explained in this paper. 

• earS. A variant processing the ears' queue while its size is greater 
than five. When the region to triangulate is a pentagon, a specific 
method using two or three incircle tests is used . Degree three, four 
and five vertices are removed by a specific algorithms. 

• flip. The flipping-based method briefly described above. 

• diirnit- A mixed method using earS if the degree of the removed 
point is > diimit and the flip method otherwise. 



Results The code was tested on a 2,000,000 points set uniformly dis- 
tributed in a square (Figure The Delaunay triangulation of the points 
is first computed (in 104 seconds) . Next the points are removed in a ran- 
dom order. We tried the methods ear3, earS, flip and the mixed method 
for 5 < diimit < 11. Figure ^ provides the whole deletion time, the number 
of incircle predicate evaluations and the number of power computations. 

These experiments have been done on a Sun UltralO 300MHz work- 
station 256Mo main memory. The code is written in C++ and compiled 



^ Filters can be classified in static, semi-static and dynamic |BBP98l: 
static: error bound is determined at compile time, 

dynamic: error bound is determined at run time, usually with an error computation for all 
intermediate results, 

semi-static: error bound is mainly determined at compile time, with addition of very few 
computation at run time. 
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Figure 3: Delaunay triangulation of random points in a square. 



with AT-T's compiler with optimizing options. Times were obtained with 
the clock command and are given in seconds and are only for the deletion 
phase. 

As infered from the theory, the performance is optimal when dumit 
is about 9. The mixed method therefore reduces the complexity of the 
deletion of high degree vertices from 0{k^) to O(fclogfc). 
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Figure 4: Deletion time for two millions of random points in a square for different 
deletion methods. Sec paragraph "Results" for details. 
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Figure 5: For ri(fclogfc) lower bound. 



5 Lower bound 

The ear-queue method, i.e. the construction of the ears in the right order, 
has a f2(fc log k) lower bound. Consider the origin r and points Pi, < i < 
n and g^, < i < n so that angles pirqi = qirpi+i — — and distances ptr = 
1 and Qir — Xi,l < Xi < a. a is chosen so that the Delaunay triangulation 
of 5 = {r, po ■ • • Pn, go . . . Qn} links r to all other points (Figure ^). When 
deleting r, all PiQiPi+i are Delaunay ears, but finding the right order on 
this ears, which is not necessary to find the new Delaunay triangulation, 
is equivalent to sorting the Xi and thus as a r2(nlogn) lower bound. 

6 Conclusion 

We have proposed a simple method for point deletion in Delaunay trian- 
gulations. This method guarantees an O(fclogfc) complexity where k is 
the degree of the removed point while most alternatives have a quadratic 
behavior in the worst case. The implementation in two dimensions corrob- 
orates these results and shows a good behavior in practice. The algorithm 
should be even more efficient in higher dimensions due to the lack of al- 
ternative methods and the higher average degree of a Delaunay vertex. 
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Code A compiled demo version is available at 
^ttp: / /www.inria.fr / prisme /logiciels / del-hierarchy / . 
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